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Abstract. We consider chains consisting of several identical subsystems weakly coupled by various types 
of next neighbor interactions. At both ends the chain is coupled to a respective heat bath with different 
temperature modeled by a Lindblad formalism. The temperature gradient introduced by this environment 
is then treated as an external perturbation. We propose a method to evaluate the heat current and the local 
temperature profile of the resulting stationary state as well as the heat conductivity in such systems. This 
method is similar to Kubo techniques used e.g. for electrical transport but extended here to the Liouville 
space. 
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1 Introduction 

As a specific topic of non-equilibrium thermodynamics, 
heat conduction has been of central interest for a long 
time. Instead of reaching a complete equilibrium state, the 
composite system under some appropriate perturbation 
enters a local equilibrium state - small parts of the system 
approach equilibrium but not the whole system. 

Within non-eqilibrium statistical mechanics the the- 
ory of linear reponse, originally developed to account for 
electric conductivity, is a very import method to inves- 
tigate dynamical as well as static properties of materials 
<Qf SHHH) I n this context the famous Kubo- formulas |H) 
have led to a rapid developement in the theoretical under- 
standing of processes induced by an external perturbation 
of the system. However, a direct mapping of these ideas 
on pure thermal transport phenomena (perturbations due 
to thermal gradients (p)) faces serious problems: Contrary 
to the case of external perturbations by an electric field, 
thermal perturbations cannot directly be described by a 
potential term in the hamiltonian of the system. Rather, 
the thermal perturbation is introduced by heat baths with 
different temperatures coupled to the system, thus calling 
for a more detailed description than is needed for elec- 
tric transport. Nevertheless, those methods are often used, 
eventually because of their immediate success in describ- 
ing non-equilibrium processes (QI K & 0) • 

Recently, the main focus of considerations on heat con- 
duction and Fourier's law has shifted towards small (one 
dimensional) quantum systems ifTol : ITl|) . Typically, these 
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systems are chains of identical subsystems weakly coupled 
by some next neighbor interaction. Based on the Lindblad 
formalism |)12|) or on techniques of quantum master equa- 
tions, heat baths are then weakly coupled to the chain at 
both ends. The influence of these heat baths enters the 
Liouville-von-Neumann equation of the system via super 
operators acting on the density operator of the chain . 
It has been found that in such systems the appearance of 
a normal heat conduction depends on the type of the in- 
teraction between these elementary subsystems ifi!}) . Most 
quantum mechanical interactions show a normal heat con- 
duction behavior (constant non-vanishing local tempera- 
ture gradients) , whereas for special coupling types the lo- 
cal gradient within the chain vanishes (divergence of the 
conductivity, non- normal scenario). 



For a comparison of these two different methods for 
investigating heat conduction - the Kubo formulas and 
the open system approach - see ijlitl . 



Instead of solving the full Liouville-von-Neumann equa- 
tion it should be highly desirable to have a simple but 
consistent perturbation method available. However, as al- 
ready indicated, the mentioned super operators describing 
the influence of the environment cannot be written as a po- 
tential term in the Hamiltonian of the system. To consis- 
tently treat heat currents and local temperature gradients 
quantitatively also for these model systems, we propose to 
extend the Kubo technique to Liouville space. 
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2 Unperturbed System 

In the following we will work in the Liouville space of the 
considered system rather than in the respective Hilbert 
space. Thus we consider super operators acting on oper- 
ators in Hilbert space, e.g. the density operator of the 
system (for a more detailed introduction to super opera- 
tors see |l5lll6Hl7j) V The complete Liouville operator of 
the system under consideration is given here by 

C = C sys + C 1 {T 1 ) + C 2 {T 2 ) . (1) 

The first term controls the coherent evolution of the quan- 
tum system defined by the Hamiltonian H: It is defined 
by its action on the density operator p according to 

C sys p = ~[H,p]. (2) 

The system H is here a chain of N identical subunits with 
n levels each, coupled weakly by a next neighbor inter- 
action, thus living in a Liouville space of dimension n 2N . 
One could think of several concrete model systems, for 
example spin models (ri = 2), for which the Hamiltonian 
would read 

N 
N-l 

(3) 

The first term is the local part of the Hamiltonian, whereas 
the second defines the interaction between the subsystems 
(crj denote the Pauli operators of the /ith spin) . Choos- 
ing J x = J y = J z we get the Heisenberg interaction and 
for J z = 0, J x = J v an energy transfer coupling only (XY 
model). Furthermore, to avoid any bias we will often use 
a random next neighbor interaction but without disorder 
(the same random interaction between different subsys- 
tems). 

The chain is weakly coupled to two heat baths, one 
at each end of the system, given by the super operators 
£i(Xi) and £2^2), where, e.g., in case of spin chains 
A(Tr)is 

£ l(Tl )p =Wj (T 1 )(2ai 1) pa« - a^a^p) 

+W}(T 1 )(2a^pa W -pa W a^-a W a^p) . 

(4) 

Here, 0+ and a_ are the raising/lowering operators act- 
ing on the first spin, and the W's are the so called rates 
which contain the bath temperature. These are standard 
Lindblad operators l)l2t) . well known from the theory of 
open systems in quantum optics, which introduce the damp- 
ing of the environement into the Liouville-von-Neumann 
equation of the system. This procedure is the same as a 



quantum master equation fliU H^) . used to describe dis- 
sipative quantum systems, for example in l)20|) . They can 
easily be generalised to finite n > 2. £2(^2) is defined 
correspondingly. 

For the unperturbed system both bath systems are at 
the same temperature T\ = T2 = T. The whole Liouville- 
von-Neumann equation of the unperturbed system then 
reads 

■^P = C Q p, (5) 

with 

C = C sys + C 1 (T) + C 2 {T). (6) 

Since the two baths have exactly the same temperature we 
expect the system to settle in a thermal stationary state 
po- This state should support neither a heat current nor 
temperature gradients - it is a global equilibrium state 
with temperature T. 

The eigenvalues and eigenvectors of the unperturbed 
system are given by the eigen equation 

C \p j )=l j \p j ) , j =0,...,n 2JV -l , (7) 

where the ket-vectors in Liouville space have been denoted 
as |. . .). A scalar product of vectors in Liouville space can 
be defined by 

(Pilfc) • (8) 

The (unique) stationary state po is also an eigenvector of 
the system with eigenvalue zero, £o\po) — 0, whereas all 
other eigenvalues have a negative real part. This is due to 
the fact that asymptotically the system should enter the 
equilibrium state \po) regardless of which state the system 
was at the beginning. No other eigenvector is able to con- 
tribute to the equilibrium state, i.e. all other eigenvectors 
must be unstable. 

Since the Liouville operator Co is not hermitian, the 
eigenvectors do not form an orthogonal basis, i.e. 

J2\Po)(Pj\=S (9) 

3 

is in general not the unit operator in Liouville space (c.f. 
ijlfih ). But with the help of the super operator Q it is pos- 
sible to find a dual basis = G~ 1 \pj) with the property 

£lft)(/^i - 1 . (io) 
3 

Based on this complete basis, we can expand each state of 
the system in terms of the eigensystem 

n 2N -l n 2N -l 

\p)= E \p3)(p j \p) = \m+ E c i\pj)> ( n ) 

with Cj = (/P|jd). Because the system will asymptotically 
be found in the normalized stationary state |po) and we 
have to require that the state \p) is normalized for the 
whole time evolution, each eigenstate must be tracefree 
(except po)). 
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3 Perturbation and Local Equilibrium State 

The system will be perturbed now by applying a small 
temperature gradient AT. We start at time t — — oo and 
switch on the perturbation exponentially till t = 0, letting 
AT constant for all times t > 0. Thereafter the system is 
subject to this small constant external temperature gradi- 
ent. We are interested in the properties of the stationary 
local equilibrium state of the system reached in the limit 
t — > oo (this state now contains stationary currents and a 
constant temperature profile). 

The Liouville operator of the perturbation is thus given 

by 

£'(AT, t) = Cx(T+ ^ff(t)) + £ 2 (T *ff(t)) , (12) 



where 



f(t) = 0{-t)c l + 0{t) , with <9(t) = 



t < 

1 t > 



(13) 

The two environment operators are the same as before but 
with a time dependent external temperature difference. 

The time evolution of the whole system under the in- 
fluence of the perturbation is given by the Liouville-von- 
Neumann equation 



^- t P=(£o + £'(AT,t))p 



Starting from the past in a thermal equilibrium state po, 
we assume the time dependent state of the whole system 
to be 

p = po + Ap(t) . (15) 

Introducing this into i|14|) , supressing terms of higher oder 
in the perturbation and observing that dpo /dt = £ p = 
one finds the time evolution equation for Ap(t), 

^-Ap(t) - £ Ap(t) = £'(AT, t)p . (16) 
at 

With a transformation similar to the one introduced by 
Kubo for the Hubert-space 

one can integrate the differential equation, finding the for- 
mal solution for the time-dependent perturbation 

Ap(t)= f dt'e^-^r'^T.OA)- (18) 

J — oo 

Let us call this equation the Kubo-formula in Liouville 
space (see also (Q)). 

In the case of a chain of two level systems each of the 
two super operators ! = 1,2 of the bath coupling at both 
ends of the system consist of two transition processes (in 
case of finite temperatures) 



with the two rates W\ (T) = (1-T)A B and W} (T) = T\ B 
(Ab is the coupling strength of the environment, T its tem- 
perature) and £\, S\ are transition operators. According 
to l@J these transition operators read, e.g., 

£{p = (2o^paf - pbfo^ - afa^p) (20) 

(For more energy levels of the subunits in the chain one has 
to account for all possible transitions - a straight-forward 
extension.) 

With this definition the perturbed super operator (|12|) 
can be rewritten as 

/\T A 

£'(AT,t) = A(T) + £ 2 {T) + —p.f{t)S , (21) 

with £ = + £\ + £2 — £^- The first two terms just 
replicate the bath operators with the same temperature at 
both ends. Acting on the unperturbed equilibrium state 
these terms vanish and therefore (|18|1 reduces to 



Ap(t) = 



ATX T rt 



&t'e Co{t - tr) f(t')£p . (22) 



Introducing the unit operator of the Liouville space, de- 
fined in JTIJ, into 



(14) A P {t) 



ATX T 



r e find using e c °^ 1 >\pj) — e lj ^ t ^\pj) 

Ap(t) = ^-J2^\ £ \M\^ J 



(23) 



M t -*'>/(t / )d* / 



(24) 

By integrating over t' with the function f(t') as defind in 
(I13fl and taking into account that the real part of lj is 
negative (for j ^ 0), we finally get the time dependent 
perturbation of the density operator 

/S«*)=^2((l+t)(p |£|«l« 



u 



I, 



(fP\£\ Po )\p 3 )\ . (25) 



£ i (T) = w{(T)£} + wJ(T)£j 



(19) 



The first term results from the integration over the addend 
j = 0, the sum contains the rest. Rewriting the matrix 
element of the super operator £ as a scalar product defined 
by ©, we find 

(p°\£\M = (P°\£M = Tr {po(£ Po )} • (26) 

The conservation of the trace of the complete Liouville- 
von-Neumann equation (|14(l requires that each operator 
produced by an action of the perturbation operator £ on 
an arbitrary density operator has to be a trace free opera- 
tor. Since \po) is the only basis state with non-zero trace, 
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as argued before, we have to require that £ po has no com- 
ponent in \po) direction. Therefore this matrix element 
should be zero and (1251) reduces to 



Ap(t) 



E 



l-h 



u 



-)(p'\£\po)\pj) ■ 

This perturbative term includes all currents and local tem- 
perature gradients of the system under the perturbation. 
Since we are interested in a local equilibrium state - a 
stationary state with a constant current and temperature 
profile, which will be reached after a certain relaxation 
time, we consider H27J) in the limit of t — > oo finding 

■n 2N -1 

Ap=hmAp(t) = — 2^ \ Pj ). (28) 



3=1 



This is the first-order change of the density operator in- 
troduced by the perturbation. 



4 Current and Local Temperature Profile 

Now we are able to account for the local temperature pro- 
file and the expectation value of the current. The energy 
current operator can be derived from a discretized 

version of the continuity equation (formulating a current 
into and out of the site p, respectively) 



i[H, ] = J 01-1 '" 5 - J ( ^ +1) 



(29) 



As a measure for the temperature T of a single subsystem 
we use here the local energy of the system so that < 
T < 0.5 in units of the local level spacing, c.f. ifbll : l2l|) . 
This should be appropriate for weakly coupled subsystems 
within the chain (only a very small amount of energy is 
within the interaction). The operator 



AH, 



loc 



H 



loc 



loc 



(30) 



measures the local energy difference between two adjacent 
subsystems p and p + 1 (p = 1, 2, . . . , N — 1). 

Since the whole stationary density operator of the sys- 
tem is now given by p = po+Ap and since we know that po 
does not give rise to any local temperature difference and 
current, the expectation values of the operators defined 
above are determined only by Ap from Ij28(l . Therefore we 
find for the local internal temperature gradient 



loc 

_\TA L! (iP\£\p ) 

3 = 1 



grpfan+i) = Tr{AH^f +1) Ap} 

" -l 



E KJ ^Tr{AHt» +1) P3} (3D 



and the local current within the system 



ATX T 



E 

3 = 1 



(P 3 |g|/5p) 

u 



Tr{J^^+ 1 )p i } . (32) 




AT 



Fig. 1. Local conductivity k< 2 ' 3 > in a Heisenberg spin chain of 4 
spins as a function of the external perturbation AT; the solid 
line refers to the solution of the full Liouville-von-Neumann 
equation, the dashed line shows 13411 . 

The current as well as the local temperature gradient are 
thus found to depend linearly on the global temperature 
difference of the bath systems. Under stationary condi- 
tions the current must be independent of p, = J, 
so that H32fl can be rewritten as 



J 



-k'AT 



(33) 



Eigenstates and eigenvalues entering here the global con- 
ductivity k' depend only on the mean temperature of 
the unperturbed system, not on AT. Based on this k' 
as a global property of the system, including its contact 
properties to the environments, let us call l|33|) "external 
Fourier's Law" . 

Furthermore, combining (|31|l and l|32|) . we can define 
a local conductivity within the system 



K 



j(fJ,,IM+l) 



j 



(34) 



implying also n^-^+i) to be independent of the external 
gradient AT. 

We can compare this result with the complete numer- 
ical solution of the Liouville-von-Neumann equation l|14l) , 
here for a Heisenberg spin chain with four spins. In Fig. [J 
we show the local conductivity of the two central spin 
systems k^ 2 ' 3 - 1 as a function of the external gradient AT. 
Indeed, we find numerically that k^ 2,3 ) according to (|34|l 
does not depend on AT (dashed line). The exact numeri- 
cal solution of the Liouville-von-Neumann equation shows 
a weak AT dependence (solid line). As expected, linear 
transport and our perturbational theory applies for not 
too large external gradients AT only. Note that the full 
range of temperatures < T < oo has been mapped here 
onto the interval [0, 0.5] (AT = 0.3 is thus already a large 
gradient). 

The theoretical predictions of the perturbation the- 
ory at hand concerning temperature gradients and cur- 
rents within the chain for all investigated systems (Heisen- 
berg, XY and random coupling model) are in very good 
accordance with the numerical solution of the complete 
Liouville-von-Neumann equation of the system (compare 

EI). 

"Normal" heat conduction (Fourier's Law) is associ- 
ated with a constant but non-zero local temperature gra- 
dient and thus a finite conductivity = n indepen- 
dent of site p everywhere in the system (see QE])). But 



Mathias Michel et al.: Heat Conductivity in Small Quantum Systems: Kubo Formula in Liouville Space 



5 



note that this internal conductivity k is not necessarily 
identical with the global k' since the latter depends addi- 
tionally on the bath contacts. In the limit of long chains 
the conductivity should be independent of the contacts, 
yielding the correct scaling behavior within the chain (see 

The majority of coupling types within the chain, the 
Heisenberg coupling and the random next neighbor in- 
teraction, indeed show this normal behavior in the weak 
coupling limit |T3: l2lh. But this "normal" transport type 
does not always show up: A spin-spin interaction consist- 
ing of an energy transfer coupling only (XY model) leads 
to a vanishing temperature gradient off the contact re- 
gions. This vanishing gradient implies a divergent conduc- 
tivity within the chain Nevertheless the current 
remains finite because of the resistance at the contacts, 
therefore the global conductivity k defined in l|33[l also re- 
mains finite for this special coupling type. Therefore we 
could state that the "external Fourier's Law" is valid even 
if Fourier's Law proper does not apply. These results can 
be compared with the numerical results of the full solu- 
tion of the Liouville-von-Neumann equation and we find 
perfect agreement. 

5 Conclusion 

We have considered heat conduction in small quantum sys- 
tems built up from identical subsystems weakly coupled by 
some next neighbor interaction. By a perturbation theory 
similar to that introduced by Kubo but extended to the 
full Liouville space of the system, we have been able to de- 
rive a quantitative equation for the temperature gradient 
and heat currents within such systems. These equations 
depend only on properties of the unperturbed system and 
linearly on the strength of the perturbation AT. 

The most remarkable point of the equation for the heat 
current and the temperature profile is the fact that the 
global temperature difference of the external bath systems 
shows up only as a parameter. This is not only a numeri- 
cal advantage, since a diagonalization for different global 
gradients is no longer necessary, but also an interesting 
physical fact: The heat conductivity is independent of the 
external gradient AT, and the "external Fourier's Law" is 
always fulfilled, even if the internal gradient of the system 
is not constant, as long as the perturbation theory applies. 

Our approach does not have the problem of introduc- 
ing a potential term into the Hamiltonian of the system, 
like in standard Kubo formulas for heat conduction. The 
bath systems, modeled by a Linblad formalism, directly 
define the perturbation in Liouville space. Like in standard 
perturbation theory in Hilbert space, the first order cor- 
rection to the stationary state of the system is expressed 
in terms of transition matrix elements of the perturbation 
operator and the eigenstates and eigenvalues of the unper- 
turbed system. Only the non-orthogonality of the eigen- 
system of the unperturbed system needs a more careful 
treatment, formally the equations are very similar. 

In the future we intend to investigate further aspects of 
the derived formalism, hoping to clarify the question of the 



different transport behavior (non-vanishing and vanishing 
local gradients) under different coupling types. 
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